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Abstract. An approximately globally convergent numerical method for a 3d Coefficient Inverse Problem for a 
hyperbolic equation with backscattering data is presented. A new approximate mathematical model is presented. 
An approximation is used only on the first iteration and amounts to the truncation of a certain asymptotic series. A 
significantly new element of the convergence analysis is that the so-called "tail functions" are estimated. Numerical 
results in 2d and 3d cases are presented, including the one for a quite heterogeneous medium. 

1. Introduction. In this paper we work with a Multidimensional Coefficient Inverse Problem 
(MCIP) for a hyperbolic PDE with the data resulting from a single measurement event. This means 
that the data are generated by either a single location of the point source or by a single direction 
of the incident plane wave. These MCIPs are non-overdetermined ones. For example, in military 
applications the best way is to collect backscattering data resulting from a single measurement event. 
This is because an installation of each new source means a life treating risk on the battlefield. 

Even though MCIPs have been studied by many researchers since 1960-ies, the topic of reliable 
numerical methods for them is still in its infancy. This is because of enormous challenges one in- 
evitably faces when trying to study this topic. Those challenges are caused by two factors combined: 
nonlinearity and ill-posedness of MCIPs. In the case of single measurement the third complicating 
factor is the minimal amount of available information. It is well known that conventional least 
squares Tikhonov functionals for MCIPs suffer from the phenomenon of multiple local minima and 
ravines. Hence, to minimize such a functional, one should apply a locally convergent numerical 
method, such as, e.g. Newton-like or gradient-like method. Convergence of such algorithms can be 
guaranteed only if the starting point of iterations is located in a sufficiently small neighborhood of 
the exact solution. However, the case when a good approximation about the solution is known in 
advance is rare in real applications. 

In a series of recent publications [3l[6l[2l|8l|9l[T0l[ni[25l[26l[28l[29l[30] the authors have 
used some properties of underlying PDE operators instead of least squares functionals. A very 
important feature of our numerical method is that it does not require any knowledge of neither 
the medium inside of the domain of interest nor of any point in a small neighborhood of the true 
solution. In all these publications convergence analysis was confirmed by numerical examples. 
Both computationally simulated and experimental data were considered. In particular, the most 
challenging case of blind real data (i.e. when the solution is unknown in advance) was successfully 
handled in [521 IMl [3D] as well as in chapter 5 and section 6.9 of the book [6J. 

For the first time, the following two goals were simultaneously achieved for MCIPs for a hyper- 
bolic PDE with single measurement data: 

Goal 1. The development of such a numerical method, which would have a rigorous guarantee 
of obtaining at least one point in a small neighborhood of the exact solution without any advanced 
knowledge of that neighborhood. 
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Goal 2. This numerical method should have a good performance on computationally simulated 
data. In addition, if experimental data are available, then this method should also demonstrate a 
good performance on these data. 

It is important to achieve both these goals simultaneously rather than just only one of them. 
Because of the above mentioned difficulties, one inevitably faces a tough dilemma in an attempt 
to achieve both Goals 1 and 2: either (1) ignore these goals, or (2) still try to achieve both of 
them. Because of this dilemma, it is natural to have the rigorous guarantee of Goal 1 within the 
framework of a reasonable approximate mathematical model (see subsection 3.5). Since convergence 
is guaranteed in the framework of that model, then we call our numerical method approximately 
globally convergent. This model is verified via a six-step procedure described in section 2. 

We are unaware about other numerical methods for MCIPs which would: (a) simultaneously 
achieve Goals 1 and 2 and, at the same time, (b) would not rely on some reasonable approximations, 
which cannot be rigorously justified. 

Compared with [31[71lllinilinilIIl[551[2nillHll2Sl[3ni, there are three main elements of this 
paper: (1) We propose a new and more convenient than before approximate mathematical model, 
(2) This model leads to a new convergence analysis, and (3) We test this model on computationally 
simulated backscattering data, both in 2d and 3d. We point out that we use the approximation of 
our model only on the first iteration of our method, see the first Remark 3.2 in subsection 3.5. 

In the majority of the above cited publications we have considered the case when the data are 
given at the entire boundary, i.e. the case of complete data collection. In the analytical study 
of this paper we also work with the case of complete data collection. However, in our numerical 
studies of section 6 we assume that only the backscattering Dirichlet data are given. Next, we use 
a numerical observation to assign the Dirichlet boundary condition at the rest of the boundary. 
Although this condition is an approximate one, numerical results show a good performance. 

The case of the backscattering data was also considered in |28j and in sections 6.1-6.7 of the 
book [6] . The 1-d case of blind experimental backscattering data was considered in [29l [30] and in 
section 6.9 of the book [6J. However, in all references cited in this paragraph both Dirichlet and 
Neumann boundary conditions were known at the backscattering part of the boundary. This led to 
the Quasi-Reversibility Method. 

In our first publications about this method the so-called "tail functions" (subsection 3.2) were 
not estimated ^iSj. Unlike this, in the new approximate model of the current paper tail functions 
are estimated on each iteration. Let C'''^" be Holder spaces, where A: > is an integer and 
a G (0, 1) . Estimating (7*^+" norms of tail functions requires all results of section 4, and this is 
the most difficult part of our convergence analysis. Indeed, we estimate functions associated with 
the fundamental solution of a certain elliptic PDE, which is valid in the entire space R^. However, 
the classical theory provides such estimates only in bounded domains [31J. Although lemmata of 
section 4 and the main theorems 4.2 and 5.1 were published in the book we believe that it is 
worthy to publish their complete proofs here as well. This is because journal publications are often 
better available and for wider audiences of readers than books. 

We search for the spatially distributed dielectric constant. We refer to [36] for a different 
numerical method for an MCIP of calculating the dielectric constant. Another non-local numerical 
method for a 2-d MCIP for a hyperbolic PDE was developed in \1M [2U] . It is based on a 2-d 
analog of the Gel'fand-Levitan-Krein equation. A different version of our approximately globally 
convergent numerical method was developed in parallel with the above publications for the case of 
a 2d MCIP for an elliptic PDE with the running source, see [27] and references cited there. The 
asymptotic behavior of the tail function in this case is radically different from ours, which led to 
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a different approximation of tail functions. This problem has an application in medical optical 
imaging of brains, see, e.g. |37) for imaging from an experimental data set for a phantom medium. 
A theory of a non-local reconstruction technique for an MCIP for an elliptic PDE with the data 
given in the form of the scattering amplitude was developed in [32^ [33]. We refer to [2] for a 
numerical implementation of this theory. Note that in numerical studies of non-local reconstruction 
techniques in [2J [Tni HO] some reasonable approximations were used, which cannot be rigorously 
justified. This is similar with our approximate global convergence concept. 

In section 2 we present the notion of the approximate global convergence. In section 3 we 
describe our algorithm. In particular, we present our approximate mathematical model in subsection 
3.5. In section 4 we prove some estimates for the function which is the Laplace transform of the 
solution of the originating hyperbolic PDE. These estimates are used in section 5 then, where we 
prove the approximate global convergence theorem 5.1, which is the central analytical result of this 
paper. In section 6 we present results of our numerical experiments. Summary is given in section 
7. 

2. Approximate Global Convergence. To verify our approximate mathematical models, 
we use the six step procedure: 

Step 1. A reasonable approximate mathematical model is proposed. The accuracy of this 
model cannot be rigorously estimated. 

Step 2. A numerical method is developed, which works within the framework of this model. 

Step 3. A theorem is proven, which guarantees that, within the framework of this model, the 
numerical method of Step 2 indeed delivers a point in a sufficiently small neighborhood of the exact 
solution, provided that the following natural condition is in place: the error, both in the data and 
in some "secondary" additional approximations, is sufficiently small. 

Step 4. The numerical method of Step 2 is tested on computationally simulated data. 

Step 5 (optional) . The numerical method of Step 2 is tested on experimental data. To have a 
truly unbiased case, blind data are preferable. This step is optional because it is usually not easy 
to actually get experimental data. 

Step 6. Finally, if results of Step 4 and (optionally) Step 5 are good ones, then Goals 1,2 are 
simultaneously achieved, and that approximate mathematical model is proclaimed as a valid one. 

It is sufficient to achieve that small neighborhood of the exact solution after a finite (rather 
than infinite) number of iterations. Next, because of approximations in the mathematical model, 
the resulting solution can be refined via a locally convergent numerical method. We have chosen 
the Adaptive Finite Element Method (adaptivity) for the latter, see [H [9l [TOl [11] and chapter 4 of 
[5]. The algorithms of our previous publications were successfully verified on two types of blind ex- 
perimental data, see [Sl[151[5ni[3n]- However, since the authors do not posses a proper experimental 
data for the algorithm of this paper, it is verified here only on computationally simulated data. 

The common perception of the term "global convergence" is that one can choose almost any 
point as the starting point for iterations, and still the process would converge to the correct solution. 
Actually, however, it is sufficient to start from such a reasonable point, which would not contain 
any information about a small neighborhood of the exact solution. In addition, it is not necessary 
to converge to the solution. In fact, it would be sufficient to reach at least one this is going along 
well with the theory of Ill-Posed problems, see, e.g. Theorem 4.6 of li4j and pages 156, 157 of |17| . 
Therefore, we come up with Definition 2.1. 

Consider a nonlinear ill-posed problem P. Suppose that this problem has a unique solution 
X* e -B for the noiseless data y*, where B is a Banach space with the norm ||-||^ . We call x* 
"exact solution" or "correct solution". Suppose that a certain approximate mathematical model M 
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is proposed to solve the problem P numerically. Assume that, within the framework of the model 
M, this problem has unique exact solution x\.j and let — x* . 

Definition 2.1 (approximate global convergence). Consider an iterative numerical method 
for solving the problem P. Suppose that this method produces a sequence of points {xn\n=i 
where the integer N e [l,oo). Furthermore, assume that this sequence is produced without any 
a priori knowledge of a sufficiently small neighborhood of x* . Let a sufficiently small number 
£ G (0,1). We call this numerical method approximately globally convergent of the level e, or 
shortly globally convergent^ if, within the framework of the approximate model M, a theorem is 
proven, which guarantees that there exists numbers Ni,N2 G [1, A^] , A^i < N2 such that 

\\xn~x*\\g<e,ne [Ni,N2]. 

Suppose that iterations are stopped at a certain number k G [Ni , N2] ■ Then the point Xk is denoted 
as Xk ■= Xgiob and is called "the approximate solution resulting from this method". 

The validity of the approximate mathematical model M of Definition 2.1 should be verified via 
above Steps 4,5. Note that the assumption of the existence of the exact x* for the noiseless data 
y* is one of the key principles of the theory of 111- Posed problems [6j [39] . 

3. The Approximately Globally Convergent Method. This method was described in 
our above cited publications. However, since we need to prove a new convergence theorem here, 
then we need to use some formulas of this method in the proof. Hence, we outline it here while still 
omitting many details for brevity. 

3.1. Statements of forward and inverse problems. Let f2 C M'^ be a convex bounded 
domain with the boundary dfl G C^. Denote \f\f^^^ = ||/|lcfc+=(o) ' ^/ ^ (^) ■ Let d = 

const. > 2. We assume that the coefficient c{x) satisfies the following conditions 

c{x)e[l,d], c (x) = 1 for X G R^\17, c G C" (r3) . (3.1) 

We assume a priori knowledge of the constant d, which amounts to the knowledge of the correctness 
set in the theory of 111- Posed problems [U [HI HZl |31] . However, we do not assume that the number 
d — 1 is small, i.e. we do not impose smallness assumptions on the unknown coefficient c{x). 
Consider the Cauchy problem for the hyperbolic equation 

c(a;)utt = AuinR^ X (0,00), (3.2) 
u {x, 0) = 0, ut {x, 0) — S {x — xq) . (3-3) 

Equation p.2p governs, e.g. propagation of acoustic and electromagnetic waves. In the acoustical 
case c{x) = b~'^{x), where b{x) is the sound speed. In the 2-D case of EM waves propagation, the 
dimensionless coefficient is c(x) = er{x), where er{x) is the spatially distributed dielectric constant 
of the medium. In the latter case the assumption c (x) = 1 for a; G R'^Xfi in p.ip means that 
we have air outside the medium of interest fl. And the assumption c (x) > 1 reflects the fact 
that the dielectric constants of almost all materials exceed the one of the air. Equation p.2p was 
successfully used in [SJ [TOl to work with experimental data, which are obviously in 3-d. The 
latter was recently explained in [12], where the Maxwell's system was solved in time domain. It 
was shown in Test 4 of [12] that the component of the electric held E {x,t) = {Ei, E2, E3) {x,t) , 
which was originally initialized, strongly dominates two other components. 
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We now formulate the CIP for the case when the data are given at the entire boundary dfl of 
the domain ft. We show in section 6 how we reduce the problem with backscattering data to this 
one. 

Coefficient Inverse Problem (CIP). Assume that the coefficient c{x) of equation Ii3.2\) 
satisfies condition iS. 1\) and is unknown in the domain fi. Determine the function c {x) for x fl, 
assuming that the following function g {x, t) is known for a single source position xq ^ 

u {x, t)=g [x, t) , V {x, t) ednx (0, oo) . (3.4) 

The function g (x, t) models time dependent measurements of the wave field at the boundary of 
the domain of interest. Practical measurements are calculated at a number of detectors, of course. 
In this case the function g {x, t) can be obtained via one of standard interpolation procedures. The 
assumption of the infinite time interval in (j3.4p is not a restrictive one, because we work with the 
Laplace transform of the function u (x, t) and the kernel of this transform decays rapidly as t — (X). 
Hence, the integral over the interval (T, oo) is actually discounted in practical computations. Thus, 
when generating the data for our CIP, we compute the forward problem for t G (0, T) , where T > is 
a finite number. Another argument here is that in our work with experimental data [51 1101 [^I291I30| 
we have actually used only a small portion of these data after the data pre-processing procedure. 

Global uniqueness theorems for MCIPs with single measurement data are currently known only 
under the assumption that at least one of initial conditions does not equal zero in the entire domain 
f2, which is not our case. All these theorems were proven by the method, which was proposed in 1981 
by Bukhgeim and Klibanov in |13[[H1[^ : also see, e.g. [151 [531 [53] for some follow up publications 
of these authors and references in books [SI [23] for publications of many other researchers about 
this method. This method is based on Carleman estimates. Actually the idea of our approximately 
globally convergent method of working with an integral differential equation, which does not contain 
the unknown coefficient c [x) , has roots in the method of [13l [lU [151 EH l22l l23l l24j . There are also 
some uniqueness theorems for MCIPs with single measurement data for the case when the unknown 
coefficient has the form const. + a (x) , where ||a|| << 1, where ||-|| is a certain norm, see, e.g. [35] . 
Our theory below does not rely on any smallness assumptions imposed on c (x). Although we image 
inclusions of small geometrical sizes in computations, the inclusions/background contrasts are not 
small. Thus, we have no choice but to assume everywhere below that uniqueness theorem is valid 
for our CIP. 

3.2. Integral differential equation. Consider the Laplace transform of the function u, 

oo 

w{x,s)— J u{x,t)e^'^''dt, ioT s > s — const. > 0. (3.5) 



We assume that the number s is sufficiently large, so that the integral (|3.5p converges absolutely 
and the same is valid for the derivatives D'^u, fc = 0, 1, 2. We call the parameter s pseudo frequency. 
It follows from (|3.2p . p.3p and p.5p that the function w is the solution of the following problem 

Aw - s^c (x) w = -6 {x - xo) , X e R^, (3.6) 
lim w (x, s) = 0, (3.7) 

|a:|— ^cxD 

see Theorem 4.1 below about p.7p . Since xq ^ fl, then Theorem 4.1 also implies that the function 
w G C^'^" (rj) . Suppose that geodesic lines generated by the function c (x) are regular and c (x) 
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is sufficiently smootli. Let t {x,xq) be tlie lengtli of the geodesic line connecting points x and xq. 
Then Theorem 4.1 of |35) implies that the following asymptotic behavior of the function w{x, s) at 
s — > oo takes place [Zl [6] 



\D':wix,s)\ 



2+a 



jjk [ exp [-ST {X, Xq)] 









} 








2+a 





, s — > cx), fc = 0, 1, 



(3.8) 



where f{x,xo) is a certain function and f{x,xo) ^ for x € fJ. It is unclear how to effectively 
verify the regularity of geodesic lines for generic functions c{x). Therefore, we assume below the 
asymptotic behavior (|3.8[) without linking it to the regularity of geodesic lines. We verify the 
asymptotics p.8p computationally, see |7| and page 173 of [5]. 

It follows from Theorem 4.1 (below) that the function w^x^s) > 0. Denote 



V (a;, s) 



Inw (x, s) 



Assuming that (|3.8p holds, we obtain 

1^ •S)l2+a = O is-^-^) , S ^ C30, fc = 0, 1. 

Keeping in mind that the source xq ^ we obtain 

Aw + (Vw)^ = c(a;), a; G 1^. 
Differentiate both sides of (I3.1ip with respect to s and let q (a;, s) — dgV {x, s) . Hence, 



V {x, s) = — J q (x, T)dT + V (x, s) , 

s 

T./ / -N \nw{x,s) 

V (X, S) = V [X, S) = ^2 ■ 



(3.9) 

(3.10) 
(3.11) 

(3.12) 
(3.13) 



Here the truncation pseudo frequency s > s is a large number. We call V {x, s) the "tail function", 
and this function is unknown. By p.lOp 



\Vix,s)\,^^ = O(s-i) , \chV ix,s)\,^^ = O (-2) ,s 



(3.14) 



The number s is the main regularization parameter of our numerical method. In the computational 
practice s is chosen in numerical experiments. 

Thus, we obtain from p. lip . (|3.12p the following nonlinear integral differential equation 



Aq - 2s^Vq j Vq {x, r) dr + 2s J Vq {x, r) dr 

s Ls _ 

s 

+ 2s^VqVV - isVV J Wq {x, t) dr + 2s {Vv f = 0, x G fl. 



(3.15) 
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By p.4p the following Dirichlet boundary condition is given for the function q 

q {x, s) — Tp [x, s) , V {x, s) e dQ X [s,s] , (3.16) 

where ip {x, s) = s^^ds In Lp — 2s^^ In </? and Lp [x, s) is the Laplace transform p.5|) of the function 
g{x,t) in (13. 4p . Equation (|3.15p has two unknown functions q and V. Therefore, to approximate 
both these functions, we approximate the function q in "inner" iterations and the function V is 
approximated in "outer" iterations, see Remark 3.1 in subsection 3.4. 

Suppose for a moment that functions q and V are approximated in ^l together with their 
derivatives g, D^V, |/3| < 2. Then the corresponding approximation for the target coefficient can 
be found via p.lip as 

c(x) = A^; + s^(V^;)^xe 1}, (3.17) 

where the function v is approximated via (|3.12p . Although any value of the pseudo frequency 
s G [s,s] can be used in p.l7|) . we have found in our numerical experiments that the best value is 
s :— s. 

3.3. Discretization with respect to s. We assume that q {x, s) is a piecewise constant 
function with respect s. Hence, we assume that there exists a partition 

s~S]\f< sn-1 < ■■■ < si < So — s, Si_i — Si — h (3.18) 

of the interval [s, s] with a sufficiently small grid step size h such that 

q{x,s) = qn{x) for s e (s„,s„_i],qo = 0. (3.19) 

We approximate the boundary condition (13.101) as a piecewise constant function, (?„ (x) = (x) , x € 
dfl, where t/i^ (x) is the average of the function ■0 {x, s) over the interval (s„, s„-i) . Next, a cer- 
tain system of elliptic equations for functions (x) is derived from p.lSp using the s— dependent 
so-called "Carleman Weight Function" exp [A (s — s„_i)] , s e (s„, Sn-i) , where A >> 1 is a certain 
parameter of ones choice. Certain numbers Ai^„, A2,n, /i,n, /o,which can be analytically calculated, 
are involved in that system, and the following estimates hold 

|Ai,«| + |^2,n| <8s2, (3.20) 
^ ifA/.>l, (3.21) 

Because of (|3.2ip we choose in our computations the parameter A >> 1 so large that we can ignore 
the nonlinear term (V^n)^ in that system. Thus, we set below 

2^ (Vqnf := 0. (3.22) 

-'0 

Our algorithm reconstructs iterative approximations Cn,i (x) G C" (flj of the function c (x) only 
inside the domain fl. To work with our algorithm, we should extend each function c„_i (x) outside 
of the domain fi. To do this, choose a smaller subdomain fi' C fi, dH.' n dfl = 0. Let the function 
X (x) be such that 

1 in fl', 

xec' (r3) , x{x) = { e [0, 1] in (3.23) 

outside of Q. 
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The existence of such functions x i^) is well known from the Real Analysis course. Let the number 
I > d. Consider the set of functions Q (d, I) C C" (O) defined as 

Q{d,l)^{ceC^ {TT):ce[i,d],\cl<l}. 

We assume in our algorithm that all functions c„ ^ G Q (d, /) . Consider the function c„ ^ (x) , 

Cn,^ (x) (1 - X (a;)) + X (a;) Cn,^ {x) , Vx e (3.24) 
Then ((X^ and ([5:^ imply that 

c„,, e ,c„,, e C" (R3) ,c(2:) = 1 for x e M'X^^- 

3.4. The Algorithm. We now describe our algorithm for approximating functions g„ and V. 
Following (jXT^ . (jXTT)) and (jXTOl) . denote 

n-l 

I'll,! (a;) = ~hqn,i (x) ~ qj (x) + Vn,i (x) ,x e il, (3.25) 

3=0 

c„_i (x) = Au„^i + (Vu„^i)^j (x) , a; € f], i = 1, ...,m, (3.26) 

where functions qj, qn,i, Vn^i are defined in this subsection below and m is the number of iterations 
with respect to tails for each given n > 1. The number m is chosen in numerical experiments. Here 
Vn^i {x) is a certain approximation for the tail function. Let Vi.i {x) be the first guess for the tail 
function, which is described in subsection 3.5. Hence, to start our iterative process, we set 

go := 0. (3.27) 

Step n^, j e [l,m] ,n > 1, see p.27p for go- For each n we iterate with respect to the tails. 
Suppose that functions qj,Vn,i (a;,s) G C^+" (il) ,j G [0,n — 1] are constructed. Then we solve the 
following Dirichlet boundary value problem for the function 

ri-l 

Ag„,i - Ain I h Vqj ) • V(?„,,; + Ai,yqn,i ■ ^VnA = 



J=0 



/«-! \^ / „_i \ (3.28) 

|90= (a;) , 

Because of p. 221) . the nonlinear term with (Vf/n.^)^ is ignored in (I3.28p . Having the function 
we reconstruct the next approximation c„_i € C"(ri) for the target coefficient using p.25p . (|3.26p . 
Next, we construct the function c„,i € C"(M3) via ([3?24|) . Next, we calculate the solution u„^i (x, t) 
of the forward problem p.2p . (|3.3p with c (x) := c„.i (x) . Next, we calculate the Laplace transform 
Wn^i (x,s) (|3.5p of the function (x,i) at s s and update the tail function using p.l3p . 

Inu-^.i (x,s) ^ r T/„,,4.i (x) if i e [1, m - 1] , , - 

32 1 (x) if i = m and n G [1, TV - 1] . 
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We set 

qn qn.n e C2+" (H) , c„ := c„,™ e (H) . (3.30) 

If i = TO and n = N , then we stop. In fact, we can stop the iterative process not only at n := iV 
but ai n := N E [li^) a-s weh. The stopping rule is chosen in numerical experiments, see section 
6. 

Remark 3.1. It is clear from (|3.28p - (|3.30p that functions qn.i are updated via inner iterations 
while "being inside" the domain Q. only. But to update tail functions Vn,i, we "go outside of f2", 
thus, using outer iterations. 

3.5. The new approximate mathematical model and the first guess Vi i [x) for the 

tail. Following the Tikhonov concept [HI 123, we assume that there exists unique exact solution 
c* {x) of our CIP with noiseless data g* {x,t) (|3.4p . We assume that 

c* G C" (R^) ,c* (x) e 1] in R^ c* (x) = 1 for x e W^\n', \c*\^ < I - I. (3.31) 

For each function c G Q {d,l) denote Wc{x,s) the unique solution of the problem p.6l) . p.7p with 
c :—c satisfying conditions (|4.2p - (|4.5[) (Theorem 4.1), where the function c is defined in p.24p . Let 
w* {x,s) be the solution of the problem (|3.6p . p.7p with c :— c* satisfying conditions (|4.2p . (|4.3II . 
Then ()4.4p and (|4.5p are also valid for w* {x,s) (see Theorem 4.1 in subsection 4.1). Using (|3.13p . 
we define tails V- {x, s) , V* {x, s) as 

lnw-(x,s) Inw* {x,s) ^ ^ _ 

V^(x,s) = ,V [x,s) = ^ ,Vs>s. (3.32) 

s s 

We call (x, s) the "exact tail". Assuming that the asymptotic behavior (|3.8p holds, we obtain 

F*(x,s) = ^^ + 0(^^^ ,s->oo,a;Ga (3.33) 

for a certain function p* (x) . We truncate the second term of this asymptotic behavior. Thus, our 
new approximate mathematical model consists of the following assumption. 

Assumption. There exists a function p* (x) E C^^" (fi) such that the exact tail function 
V* {x, s) has the form 

V* {x, s) Furthermore, '—^ = Vs > s. (3.34) 

s s s 

Since g* (a;, s) = dgV* {x, s) for s > s, we derive from (|3.34p that 

q* ix,s) ^ (3.35) 
s 

Set in p.l5p s = s. Then, using p.34p and p.35p . we obtain the following approximate Dirichlet 
boundary value problem for the function p* (x) 

Ap* = in O, p* e C2+" (n) , (3.36) 
P*\dn = -s'^i^* {x,s) , (3.37) 
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where 'ip*{x,s) is the exact function 'ip{x,s), which corresponds to the function g* (x,t) . The 
approximate equation (j3.36[) is valid only within the framework of the above Assumption. Although 
this equation is linear, formula (|3.17p for the reconstruction of the target coefhcient c* is nonlinear. 
Recall that by (|3.16p q (x, s) — (x, s) , V {x, s) e dfl x [s,s] . Assume that 

(x, s) e C^+" (H) , Vs e [s, s] . (3.38) 

Consider the solution p (x) of the following boundary value problem 

Ap = Qmn, pe C2+" (H) , (3.39) 

p\on = -s^^{x,s). (3.40) 

As the first guess for the tail function we take 

1^1,1 {x) := (3.41) 
s 

By the Schauder theorem there exists unique solution p of the problem (|3.39|) . (|3.40p . Furthermore, 
it follows from (|3.36p - (|3.4ip and Schauder theorem that with a number M — M (ft) > the 
following estimates hold 

|VFi,i - Vy*|,+„ < MU{x,s)-r {x,s)\\c^^<^i9n) , (3-42) 
|VFi,i|,+„<Ms||V(x,s)||c,2+.(gf,). (3.43) 

Remarks 3.2. 

1. The main approximation is the second equality p.34p . This approximation amounts to the 
truncation of the second term O (s^^) of the asymptotics p.33p . It is made only to obtain the 
estimate p.42[) on the first iteration of our method to obtain estimate (|5.36p in the proof of the 
convergence Theorem 5.1. On all follow up iterations in that proof we do not use the second equality 
(|3.34p . Rather, we use the true fact that V* (x) — Inw* {x,s) . 

2. It follows from (1?^ that, substituting (|XIT|) in (j?^ and (jX^ at n = 1 and setting 
91,1 ■= 0, we obtain a good approximation for the exact solution already on the first iteration of 
our method, as long as the error in the boundary data tp* {x,s) is small. The smallness of the error 
is a natural assumption. Theorem 5.1 guarantees that all other solutions obtained in the iterative 
process of subsection 3.4 also provide good approximations, as long as the number of iterations is 
not too large. This means that we should develop numerically a stopping criterion to stop iterations, 
see section 6. Suppose now that iterations are stopped before this stopping criterion is met, e.g. 
just on the first iteration. In this case we can apply the second stage of our two-stage numerical 
procedure [6l HI |9l [TOl [TTj . Namely, we could apply a locally convergent numerical method to refine 
the solution via taking the solution obtained on the globally convergent stage as the starting point 
of iterations. Such a method can be applied indeed, since Theorem 5.1 guarantees that the iterative 
solution, at which we have stopped, is close to the exact solution c* {x) . Numerical confirmations 
of this can be found in tests 2,3 of |H] and in tests 2,3 of section 4.16.2 of [Ij. 

We now establish uniqueness within the framework of our approximate mathematical model. 
We refer to Lemma 2.9.2 of [6] for the proof of Lemma 3.1. 

Lemma 3.1. Let the above Assumption holds. In addition, let ^3. 1 7]) holds for the function 
V* {x,s) , 

V* (x, s) = — / q* (a;, r) dr + (x, s) 
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with the tail function V* {x,s) satisfying conditions ^3.34^ , 113.36]) and 113.37]) . i.e. 

c* (x) = \Av* + s^ \Vv* 1^1 (a;, s) , (x, s) e f7 x [s, s] . 

Then there exists at most one function c* (x) . 
4. Estimates of Tails. 

4.1. Some estimates of the function w; (a;, s). First, we should justify (|3.6|) . (|3.7p . Theorem 
4.1 is a combination of Theorems 2.7.1 and 2.7.2 of [6]. Therefore we refer to [6J for the proof. 
Theorem 4.1. Assume that the coefficient c{x) of equation ^3. S]) is such that 

c e C''+" (M^) ,c(a;) e [l,d], c (a;) = 1 for a; G R^\r2. (4.1) 

Assume that the solution u {x,t) of the problem ^3.'^) . i3. 3]) is such that 

u£C^{t> Ml {x,xo,cj) ,\D'^u{x,t)\ < Ah {xo,c) ef^Wj] < 2,Va; e R^yt > Mi (x,xo,c), 

where Mi{x,xo,c) and /3 — 13(c) are positive numbers depending on listed parameters. Let the 
function w (a;, s) be the Laplace transform ^3.5\) of the function u {x, t) . Then there exists a number 
s — s{c) > (3 (c) such that for all s > s (c) the function w (a;, s) is the unique solution of the problem 
of the form 

, cxp (— s la: — a;o|) \ 1 1 n\ 

w (x, s) = — j j h w [x, s) ■— wi (x, s) + w (x, s) , (4.2) 

Att |a: — a;o| 

w £ c'=+2+" (rS^ (4 3) 

Also, the following inequalities hold 

Wd (x, s) < w (x, s) < wi (x, s) ,\fx ^ xq, (4-4) 
where the function {x, s) is the unique solution of the problem i3.6]) . 113. 7)) for the case c (x) = d, 

exp ( — s^/d \x — xq\ 



Wd{x,s) = — ■ . (4.5) 

An \x — xo\ 

Furthermore, consider the problem 113.6]} . i3. ?)] irrelevantly to the problem 113.2]} . 113.3]} while still 
assuming ^.1\ }. Then for any s > there exists unique solution w {x, s) of this problem satisfying 
conditions ^.2^ , i4.3^ . Furthermore, conditions |4.^[ ), ^.5\ } hold for this function w (x, s) . 

Lemma 4.1. Let x{x) be the function defined in 113. 23\) . functions Ci,C2 £ Q{d,l) and func- 
tions ci (x) ,ci (x) be defined as in ^3.24\ }. Then \ci — C2\^ < \x\a 1*^1 ~ ^'^\a ■ 

Proof. By (|3.24p ci (a;) — C2 (x) = x (x) (ci (a;) — C2 (a;)) . The rest of the proof follows from 

\f9L<\f\jgL,yf,g£C-{n).a 

Note that (il) C C" (fi) , and also there exists a constant C — C (il, a) > such that 

l/L<c||/||^,(^),v/eci(n). (4.6) 

Lemma 4.2. Let the source xq ^ Q.. Then there exists a constant Y = Y (fi, s, d, I, x,Xo,a) > 
depending on listed parameters such that \u>c ix,s)\^ <Y,\/c £ Q (d, I) . 
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Proof. Below in this proof Y ~ Y {il, s, d, I, x, xq, a) > denotes different constants depending 
on listed parameters. By (|4?2)) and ((43| ^Pc {x, s) e 6"^+" (H) . Denote 6 (x) = c (a;) - 1. Then 



w-i;ix,s) = wq{x,s) ~ s J wi{x-^,s)b{Owc{^,s)d^. (4.7) 
o 

By gll) and gT]) 

(x, s)| <y + y J wi{x-^,s)d^<Y, xen. (4.8) 
n 

In addition, by (|4Jl) 

Vu>c (x, s) = Vwq (x, s)-s'^ J Vwi (x - ^, s) 6 (C) u;- (^, s) d^, x e fi. (4.9) 

Hence, |Vw;c(x,s)| <Y, x efl. Hence, (|46|) . ([48| and ([49| imply that lit^ (x,s)|^ < F. □ 
Consider a bounded domain ili C M'^ such that 

n c r^i, dfi n af^i = 0, drii e c^, xq ^ Hi. (4.10) 

Lemma 4.3. Let 17' C il C Oi be above bounded domains in M"^, condition (^T7^ 6e satisfied 
and x{x) be the function in 113.23]} . Let s > 1. Then the function Wc{x,s) £ (dili) ,yc £ 
Q{d,l). Furthermore, there exists a constant B = i? (f2, fi', r2i,s, d, x, Xq, a) > 2 depending only 
on listed parameters such that 

\\w^ix,s)\\c.^on,)<B, yceQ{d,l). (4.11) 
For any two functions Ci, C2 G Q {d, I) denote w (x) = Wc-^ (x, s) — Wc^ {x, s) . Then 

\\w\\c^9n,)<B\ci~C2\^, yci,C2eQid,l). (4.12) 

Proof. Everywhere below in this paper B denotes different positive constant depending on 
above parameters. The integrand of formula (j4.7p does not have a singularity for x G f2i\r2. Hence, 
(|4Jl) implies that Wc{x,s) G {dili) . Next, (|44l|) follows from (|4J|) and Lemma 4.2. 

Denote 

c (x) = Ci (x) - C2 (x) , 61 (x) = ci (x) - 1, 62 (a;) = C2 (x) - 1. 

Hence, by p.24p ci (x) — C2 (x) = x (a;) c (x) . First, substitute in (|4.7D (foi, lUcJ • Next, substitute 
(&2,Wc2) • Next, subtract the second equation from the first one and denote w (x) = Wc-^ {x,s) — 
{x,s). We obtain 

w (x) = ~s^ fwi{x- C, s) X iOciO w-c, (^, s) - s2 / wi (x - e, s) 62 (0 55 (0 (4.13) 
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Let 

h (x) - -s^ / _ s) X {£.) c {0 w-c, , s) d^, 



-2 



a 

Using the same arguments as ones in the proof of (|4.1ip . we obtain 

ll^i|lc3(ao,) <^|21o- (4-14) 

Next, 

Aw — s'^C2W — s^x {x)c{x) uJci {x,s) , x G R'^. (4-15) 

Fm'thermore, it follows from (|4.13p that the function w (x) decays exponentially together with its 
derivatives as \x\ — > oo. Hence, multiplying ()4.15p by w, integrating over M.^ and using Lemma 4.2 
and the fact that x (x) = for x G ]R'^\f2, we obtain in a standard manner ||'5'||x,2(n) — ^ ll^lL2(f2) — 
B |c|^ . Hence, ||/2|lc3(aoi) ^ ^ l^a ■ This estimate combined with ()4.14p implies ()4.12p . □ 

4.2. Estimates of tails. Theorem 4.2. Let il' C C J^i C M"^ be above bounded domains 
with condition and let x{x) be the function in 13. 23]) . Also, let s > 1. For each function 

c £ Q {d, I) consider the function c defined in Iji3.24.\ ). Denote 

V.(x)^'^^^^^. (4.16) 

Then there exists such a constant B ~ _B (fJ, 57', r2i,s, d, xq, a) > 2 depending only on listed 
parameters such that 

|vyeli+„<s,Vceg(d,o, (4.17) 

iV^ei - VycJi+„ < S |ci - ,Vci, C2 e Q (d, . (4.18) 



Proof. By 



VV.(.) = OlMx) = - = 1,2,3. (4.19) 

SW-c{x,S) SWc(X,S) sw^(x,s) 



Since by (|44l) and ([431 



4^exp(s|2;-:Eo|) ^ , ^ 1 ^ 47r exp (sVd |a; - 

s^w- (a;, s) 



k - a;o| < .^2 ^^ < =2 - a;o| < B, 



then (|4.19p implies that in order to prove (|4.17p and (|4.18p . it is sufficient to prove that 

\w-{x,s)\^^^<B,yceQ{dJ), (4.20) 
\wc, {x,s) - Wc^ (a;,s)|2+„ <B\ci- C2U • (4.21) 
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Denote fc{x) — 'Wc{x, s) jaoj. By Lemma 4.3 



(4.22) 



Since by (jlT^ x^i ^ VLi, then and imply that (a;, s) S C2+" . On the other hand, 
the function 'w-{x,'s) solves the following Dirichlet boundary value problem in Hi 



Ai«c — s^c (x) Wc = 0, X e r^i, 

u'c I OHi /- (a;) • 

Since the C^+" (fii) solution of this problem is unique, then Schauder theorem and (|4.22p imply 
that (a^, s)l2+a — B Il/llc2+t«(ani) — which proves (|4.20p . Denote again w {x) = {x,s) — 
{x,s) . Then 

Aw — S^Ci (x) W = s'^ (ci (x) — C2 (x)) Wc^ , 

w I aoi = /ci {x) - fc2 (x) ■ 

Hence, Schauder theorem and Lemmata 4.1 and 4.3 imply (|4.2ip . □ 

5. Approximate Global Convergence of the Algorithm of Subsection 3.4. 

5.1. Exact solution. Recall that we the existence and uniqueness of the exact solution c* (x) 
of our MCIP satisfying p.3ip . The corresponding functions w* {x,s) , V* {x,s) were defined in 
subsection 3.5. Let s > s (c*) , where s(c*) is the number of Theorem 4.1. Let u* {x,t) be the 
solution of the problem (13. 2p . (13. 3p with c :— c* . Since one can differentiate infinitely many times 
with respect to s under the integral sign in p.Sp . then 



Denote 



1* (a^.s) = ^ 
OS 

Consider functions q* {x) , -ip^ {x) 



w* {x, s) e C2+" (n) X [s, s] , if s > s (c*) . 
d [In [w* {x, s)] 



, -0* {x, s) = q* {x, s) \on,s e [s, s] . 



(5.1) 
(5.2) 



9n(a;) = T / q* ix,s)ds, ipnix) = - / ip* {x,s)ds, Qoix) = 



(5.3) 



Then (|5lJ) and (|0|) imply that 

k* (x, s) - ql (x)|2+„ < C*h, ij* (x, s) - ^* (x) 



< C*h,ne [l,N],se [s„,s„_i] 



(5.4) 

Here the constant C* = C* (ll'7*llc2+a(n)xciLs s] > depends only on the C^+" (fi) x [s,s] 
norm of the function q* (x, s). Hence, we can assume that 



jnax k:i2+„ < C* 



l<n<N 



(5.5) 
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Without any loss of generality we assume that 

C* > I. 



(5.6) 



By the one of concepts of Tikhonov (see, e. g. section 1.4 of |6]) we assume that the constant C* 
is known a priori. By (|5.2p 

q: (x) - C (^) , ^ e dn. (5.7) 
Hence we obtain the following analogue of equation p.28p 



3=0 
n-1 



J=0 



(5. 



2^2,nVy* U ^ Vg* - ^2.n | VF* f + Fn (x, h, A) . 



3=0 



Here the function F„ {x,h,X) S C" (57). The term 2/i^„ (Vg*)^ //q is included in i.e., unlike 
p.22p . we do not ignore this term now, since we work now with the exact solution. Hence, by p.2ip 



Let 



max|F„(a;, /i,A,s)|„ <C*h,n€ [l,N]. 

Xh>l 



ix)^-hq:{x)^hY,qj{^) + V*{x), xen, ne[l,N]. 

3=0 



Using (|5.4p . we obtain similarly with p.26p 



c* (x) 



(x) + Fn (x) 



where the error function Fn is such that 



\Fn\ <C*h. 



(5.9) 



(5.10) 



(5.11) 



(5.12) 



We also assume that the function g{x^ t) in (|3.4p is given with an error. This naturally produces 
an error in functions -!/;„ in (|3.28p . Let cr > be a small parameter characterizing the level of the 
error in the data tp {x, s) . Because of p.38p . we assume that in p.28p functions ?/;„ (a;) G C^+" (9$!) 
and 



<C*{a + h). 



(5.13) 



c2+°(af2) 

In addition, we assume that We now reformulate the estimate of the Schauder theorem for the 
specific case we need. Consider the Dirichlet boundary value problem 

3 

Au + '^bj{x)uxj - bo{x)u ^ f (x) ,x e il, (5.14) 



3=1 



u I an^gix)^ [dn) 



(5.15) 
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Assume that the fohowing conditions are in place 

bj,bo, f eC" (H) , 60 (x) > 0, max {\bj\) < P, P ^ const. > 0. (5.16) 

jG[0,n] 

Then Schauder theorem |3T] claims that there exists unique solution u e C^^" (fl) of the bound- 
ary value problem (|5.14p . (|5.15p . and the following estimate holds with a certain constant K = 
K {Q., P) > 2, depending only on the domain and the constant P 



9\y^Han) + \fL ■ (5-17) 



5.2. Approximate global convergence theorem. Theorem 5.1. Let O' C O C fii C 

be above bounded domains with condition ^.10\ l, let x (x) be the function in US. 23\) . s > 1 and 
^3.38]} be valid. Let the function c* {x) satisfying conditions fi3.31\) be the exact solution of the 
CIP fi3.}A] - ^57^ , where constants d,l > 1 are given. Also, let condition i5.13\) holds, where a 
is level of the error in the data, in i3. 28\) functions i/j^ G C^+" (dfl) , and the constant C* = 

^* (ll9*llc2+°(r2)xci[s s] ^ 1 defined in ^5.4^ - [57d\) . Consider the algorithm of subsection 3.4 
supplied by Assumption of subsection 3.5. Let the first tail function Vi_i (x) be calculated via i3.39\) . 
^3.40^ and \3.4i\ l. In addition, assume that all functions Cn,i (x) in i3.26]) are such that 

Cn,tix) >l,x en. (5.18) 

Assume that the parameter A of the Carleman Weight Function is so large that Aft. > 1 (see i3.21\} ). 
Let I15.13\) be valid and also 

IIV- {x,s) - 1})* (x, s)||c.2+c(an) < C*a. (5.19) 
Consider the error parameter rj, 

r] = h + a. (5.20) 

Let B =^ _B (ri, s, d, Z, X, 2^0, ck) > 2 be the constant of Theorem 4-2. Consider the number Bi = 
Bi {fl,D,i,s,dJ,C* ,x,xo,a) , 

Bi = max (4B + 3C* , 245^) > 24. (5.21) 

Let 

K ^K{fl, i^Bi) > Bi (5.22) 
be the constant in 115. 17\} . Let the parameter rj be so small that 

^e(0,%), rfo^ ^^^3^^ . (5.23) 

Then 

c„,, e (H) ,cn,^ e (R3) , [n,i) e [1, A^] x [l,m] , (5.24) 

c„,j,c„,, {x) e Q {d,l) , {n,i) e [1, A^] x [l,m] . (5.25) 
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In addition, the following estimates hold for {n,i) G [1,-^] x [1,to] 

|VKmIi+„<Si, (5.26) 
|VK,. - yV*\,^^ < s3[»-i+(n-i)H+i . (5 

\qn,-q:\,^^<KBf^^-'^"'^-rj, (5.28) 
|gn,^l2+„,|g„l2+a <2C*, ne[l,N], (5.29) 

|cn..-c*|„<i3f+("-^'"1.77. (5.30) 

Define the number uj € (0, 1) as 

In (iCA^) 



2[3iVTOlnBi +ln(ifiV)] 



,c<.e(0,l). (5.31) 



Then \5.80\) becomes 



<t;" :=£€ (0,1). (5.32) 



Therefore, by 115.321) and Definition 2.1 the algorithm of subsection 3.4 possesses the approximate 
globally convergent property of the level e. 
Remarks 5.1: 

1. Since K — K{n,s'^Bi) , one can incorporate the term KN in (|5.23p in the term Bf^"^. 
However, we are not doing this for the convenience of the proof. 

2. Condition (|5.23p provides a hnkage between the level of the error rj in the data and the 
total "allowable" number of iterations Nm. The fact that the maximal number of iterations Nm is 
limited is going along well with the theory of Ill-Posed Problems. Indeed, it is well known that the 
maximal number of iterations and the error in the data are often connected with each other, see, e.g. 
pages 156 and 157 of [T^ and section 1.6 of jB]. Hence, So that the maximal number of iterations 
Nm is a regularization parameter in this case. The fact that the constant Bi depends not only 
on the domain Q but also on the domain Qi does not affect the approximate global convergence 
property. 

3. It is hard to establish a priori the upper limit for the maximal number of functions (7„.We 
have consistently observed in our numerical tests that certain numbers indicating convergence sta- 
bilize a few iterations before a certain number iV e [1,^] , i-e. at a certain n < N. Next, they 
grow steeply for n^> N. This means that the process should be stopped at a certain n = N < N. 
Usually we take N := N — 1, see, e.g. pages 178-182 and 311-314 in the book [6J. This numerical 
observation is going along well with (|5.23p . (|5.32p . 

Proof of Theorem 5.1. Estimate (|02|) follows from estimates ([5?23l) . ([O0|) and ([OT|l . 
Hence, we focus below on the proof of relations (|5.24p - (|5.30p . Denote 

VnA = Vn,i -~ V*, QnA QnA ^ Qn, Qn = Qn ~ Qn, 

By p.23p . (|3.24p and (|3.3ip c* (x) = c* (x) . Hence, we can apply Theorem 4.2 to estimate the norm 
Suppose that estimate (|5.30p holds. Then Cn,i S Q {d,l) . Indeed, using p.3ip . (|5.23p 



l+a 

and (|5.30p . we obtain 



C* +CX< |C*L + \Cn., -CX<1-1 + . ^ < L 
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Similarly Cn,i < d. The latter two estimates and (|5.18p imply that c„,i £ Q {d, I) . Next, since the 
function Cn.i E Q {d,l) , then, using (|4.17|) and (|5.21l) . we obtain (|5.26[) . Also, since the function 
Cn,i G , then the function Cn,i & [^,d]. Hence, if (|5.30p is true, then (|5.24p and (|5.25p hold. 

First, we prove (|5.28p - (|5.30p for the case (n,«) = (1, 1) . Subtracting equation (|5.8p from equa- 
tion p.28p at {n,i) — (1, 1) and also subtracting (|5.7p from the boundary condition in p.28p . we 
obtain 

Agia + Ai,iVyi,iVgi,i = -Ai.iVi^iaVgt - ^2,iVi^,i (V(/i,i + VF*) - i^i, (5.33) 
qi 1 (a;) = V'l (a;) ,x e 917. (5.34) 

We now estimate the right hand side of It follows from ((X?T|l . ([g:^ and (jiTf)) that 



By (lOl) and (l5T9ll 



l+a 



< B?7. 



(5.35) 



(5.36) 



Estimates ([OS)) . (I07|) for (n, i) (1, 1) with Bi being replaced with B follow from (gSll), (gSHl)- 
dSHI). Using (1131), (IE211), ((071) for (n, i) = (1, 1) with Bi being replaced with B, as weU as ([2201), 
dSini), (EH), (|OT|) . (|05)) and dOH), we obtain 



Ai,iVFi,iVgi* + (VFi,i + VV"*) - Fi 

< 8s^BC*r] + Ws'^B^T] + C*r] = Ss^B \ 2B + C* + 



C* 



?7 < As^BBiT] < s'^Bfrj 



Thus, 



Ai.iVt/iaVgi* + A2,iVT4,i (VVi,i + VV*) - Fi 



< s^Bfr]. 



(5.37) 



By dEHl) 



C* < 



Bi 



(5.38) 



Next, consider coefficients in the left hand side of equation (|5.33p . Using (|3.20p as well as (|5.26p 
for (n,i) = (1,1), we obtain |Ai_iVVi^i|^ < 8s^-Bi. Hence, conditions (|5.16p are satisfied with 
P = Ss'^Bi. Hence, by (|5.17p the solution of the Dirichlet boundary value problem ()5.33p . (|5.34p 
can be estimated as 



\qiM2+a<^'KBlr^ + K 



^1 



c2+°(an) 

Using (|5.13p . (|5.20p and (|5.38p . we obtain from this inequality 



,K (n,s'^Bi) = const. > 2. 



\qi,i\,+^<KBt[s 



C* 



f] < KBj 



1 

72 



77 < 2s'^KBf'q. 
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Hence, applying (|5.2ip . we obtain 

l9i.il2+a<^^iV (5-39) 

Estimate ([QS]) for (n, i) = (1, 1) follows from ([09ll . Next, using JSJ]), (|03l) and ((09| . we obtain 
(1?:^ for {n,i) = (1,1), 

|9l,ll2+a < I9ial2+a + l9ll2+a < ^^^^ + C* < 2C* . (5.40) 

Since by ((3:391) and (IXiTj) ^1,16 (II) , then (13:25)) and (lOO)) imply that U14 e 6*2+" (H) . 

Hence, by ((X^ cij e C" (H) . This, ((X^ and ((X^ imply that ci,i e C" (R^) . Hence, ([Oi)) 
is true for {n, i) — (1, 1) . 

Now we estimate the norm |ci . Subtracting (IS.lip from (|3.26p for (n, i) = (1, 1) , we obtain 

ci.i = Avi^i + 4 Vwi,i (Vvi.i + \7vl) - Fi. (5.41) 

Hence, using (|5J2l) . ([QO]) . (|08l) and (lOTj) . we obtain 

l^i.iL < |V?^i.ili+„ [l + s' (|V^;i,iL + |V«i*L)] + ^V- (5.42) 

Subtracting (lETU)) from we obtain = + Fi^. Hence, it follows from ([05)) . ([07|) 

at (n, i) = (1, 1) and from (I5.39P that 

\^vi,i\i+a < KBlif + B177 < 2Bi77. (5.43) 
By dESD, (QUI) . (lOTD . (I05D and (I05D - (I05D 

|V<|i+„<C*7Vr, + Bi<2i3i. (5.44) 
Hence, (fOSl) . ((Oil) and ((OS)) imply that 

|Vi;i,i|i+„ = |Vvi,i + Vi;i*|i+„ < |Vwi,i|i+„ + 2Bi < 2Bir] + 2Bi < Wi. 

Hence, 

l + s^dVi-ial^ + lVt-ri^) < l + 5Bis2 <6s2Bi. (5.45) 

Hence, comparing ()5.45p with (|5.42p and (|5.43p . we obtain 

|ci,iL < 12s' 5??? + C*77 < Blr^. (5.46) 

This establishes (|5.30p for (n,i) = (1, 1) . As it was proved above, (|5.30p for (n, i) = (1, 1) implies 
((5:241) and ((OSj) for (n, j) = (1, 1) . In summary, we have established ((5:24]) - ((OOj) for (n, i) = (1, 1). 

Since we have established relations (I5.24l) - (|5.30l) for (n, i) = (1,1), we can assume now that 
we have proved (|5.24p - (|5.30p for {n',i') £ [0,n] x [0,i— 1], where n > l,i > 2. We now want to 
prove (|5.24p - (|5.30p for {n',i') = {n,i) . The mathematical induction principle and formulas p.29p 
and p.30p imply that this would be sufficient for the proof of Theorem 5.1. 
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Subtracting equation (j5.8p from equation p.28p and taking into account boundary conditions 
(see (|5.7|) ). we obtain 



Ai,„V(?; - A2,nh ^ (Vg, + Vq*) + 2A2,„VK,« U I] Vq, 



J=0 



3=0 



(5.47) 



First, we estimate the difference of tails Vn^i- Since relations (|5.24p - (|5.30p are valid for {n',i') G 
[0, n] X [0, 1 - 1] , then by Theorem 4.2 

VF„,.|i+, <i3i, (5.48) 



l + Q 



l + Q 



These estimates establish (|5.26p and (|5.27p for (n',i') = {n,i) . 

We now estimate the right hand side of equation (|5.47p . First, using p.20p . (|5.5p . (|5.35p and 
(|5.48p . we obtain 



n-l 



This inequality and (|5.38p lead to 



< 8s^ (C* + 3C*Nh + 2Bi) < 8s^ (C* + 1 + 2Bi) 



3=0 



(5.50) 



Estimates (|5.28p hold for functions qj = qj — q*,j E [0,n — 1] . Hence, using (|5.23p . we obtain 



3=0 



< KBl^"'NTf < r] 



Combining this with (|5.50p . we obtain the following estimate for the term in the second raw of 

(EUD 



3=0 



hY.^l3 



3=0 



< 24s^Bi?7. 



(5.51) 
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Next, using (jX^ . (jS^]) . (|OT|l . (jO^ . and ([SlIH)) . we obtain 

Hence, using (jOl) . (|5.38p and (|5.49p . we obtain 



\Fn\ 



3 

Combining this with (|5.5ip and using (|5.2ip . we obtain the following estimate for the right hand 
side {rhs) of (fOT]) 



\rhs\^ < 21s^BiB'l^ 



< 21s2BiBf'"'+^""'''"l+' 



24 
2IS1 



1 + J_)^^22s2BiBf-i+("-i)™l+i 



Thus, 



(5.52) 



We now estimate coefficients which are multiplied by Vq„^i in the left hand side of (I5.47p . We 
use p.20p . ()5.2ip . (|5.22p . ()5.23p and (|5.38p . By the assumption of the mathematical induction 
method we have that inequalities (|5.29|) are valid for functions qj with j G [0, n — 1] . First, 



j=0 



Next, using IfJ^ and we obtain 



6s^ 1 
Bf^ < J^- 



(5.53) 



(5.54) 



Hence, it follows from (|5.53p and (|5.54p that the Dirichlet boundary value problem (|5.47p satisfies 
conditions (|5l4l) - ([5T7l) with P ^ 9s^Bi,K ^ K {n,s'^Bi) > 2. Hence, using ([ETO)) . ([ET7|) . 
and ([53^ . we obtain 



I2+Q 



22s2SiSf'"^+^""^^'"l+^ + C* 



Since by (|5.2ip Bi > 24s^, then the last estimate leads to 



|gn,4+„ < KB, 



3[i+(n-l)m] 
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a) G — Gfem U Gfdm b) Gfdm c) Gfem — ^ 



Fig. 6.1. a) Geometry of the hybrid mesh. This is a combination of the quadrilateral mesh in the subdomain 
Gfdm b), where we apply FDM, and the finite element mesh in the inner domain Gfem = ^ c), where we use 
FEM. The solution of the inverse problem is computed in Gfem = ^■ 

which is dOS). Next, we prove We use (H^]), JESl), ^f^^ and 

Estimate now the norm |c„.i|^ . We obtain similarly with (I5.42p 

|2n,4 < [1+^2 {\VvnA^ + |V<|J] + l^r;. (5.55) 

We have 

n-l 

i=o 

Hence, by jS^Ql), (lOT]) . ([OSl . ([ST^fll and ([QS]) 

|V«„,|i+„ < if7Vij3[^+(«-l)H^2 ^ ^3[,-l + („-l)™] + l^ ^ 25^3[.-l+(„-l)™] + l ^_ ^5 

Next, using (|3.25l) and (|5.44p . we obtain similarly with (|5.45p 

l + s'(|Vt.„,4 + |V<U) <6s2Bi. 
Combining this with (|5.55p and (|5.56p and using (|5.2ip . we obtain 

< 52^3[j-l + («-l)m] + l ^ ^ ^ ^3[i+(n-l)m] ^ ^ 

Thus, |c„,,;|^ < 53[j+(n-i)m]^^ .^^i^jgi^ provcs ([OO)) . Tlius, relations ([E^ - ljOn]) are vahd = 
(n,j).n 

6. Numerical Studies. In this section we conduct some numerical experiments in both 2d 
and 3d cases. In the 2d case we use specific ranges of parameters for a simplified mathematical 
model of imaging of antipersonnel land mines, see [53] and sections 6.8.2 and 6.8.3 of [6J for this 
model. In the 3d case we model imaging of explosives hidden on belts worn by humans. We point 
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out that in both cases our mathematical models are certainly simplified ones and further studies 
are necessary to see how they reflect the reality. 

It is well known that there are always some discrepancies between the theories and numerical 
implementations of complicated numerical methods. We now list two discrepancies for our case. 
First, the above theory was developed for the case of the point source, because of a convenience of the 
analysis. In computations, however, we work with the case of an incident plane wave with the single 
direction of incidence. This is because it is better to operate with a plane wave computationally. 
Also, in the case when the point source is far from the domain of interest, it can be approximately 
treated as a plane wave. The above theory can be extended to the case of a plane wave after 
a purely technical additional effort. Second, to decrease the complexity of our computations, we 
replace p.23p and (|3.24p with the following simplified formula 

_ ( \ _ j c„,i (x) if Cn,i (x) > 1 and a; e O, , . 

'''''' - \ 1 if either c„,, (x) < 1 or x ^ U. ^ ' 

When reconstructing functions Cn,i (x) , we use a weak formulation of (|3.26|) via finite elements, 
see pages 184, 185 of [6] for this formulation. Comparison of Figures 3.12 and 3.13 of [6] (pages 
182, 183) shows that this formulation provides significantly more accurate results than the strong 
formulation (|3.26p . 

We now describe our stopping criterion used in computations of sections 16. 1[ 16.21 We stop 
computing functions c„_i on every pseudo-frequency interval [s„,s„_i) when 

either 7V„ > N^-i or 7V„ < 77, (6.2) 

where 

AT I |Cn,i ~ C„^i_i I o\ 
Nn = (6.3) 

l|Cn,i||L2(n) 

Here, i is the number of iterations with respect to the tail on every pseudo-frequency interval 
[sn, Sn-i)- Recall that we define by m the number when iterations with respect to the tail are 
stopped. 

To generate data for the CIP, we solve the forward problem for equation p.2p . Since it is 
impossible to numerically solve this problem in the entire space {n = 2,3), we solve it in a 
rectangle in 2-d and in a rectangular prism in 3-d, just as in [B]. We denote this each of these 
domains G. Thus, G is our computational domain in which we compute the forward problem, and 
it replaces M" (n = 2,3) . We impose the first order absorbing boundary condition [TH] on one part 
of the boundary dG and zero Neumann boundary condition on another part of dG. In all cases the 
domain of interest il C G, dft n dG = 0, see for details below. 

6.1. Our mathematical model of imaging of plastic antipersonnel land mines: 2d 
study. The first main simplification of our model is that we consider the 2d case instead of 3d, 
although a 3d numerical test is also presented below. Second, we ignore the air/ground interface, 
assuming that the governing PDF is valid on the entire 2d plane. Results of studies of experimental 
data in |291 [30] as well as of section 6.9 of [6J indicate that the influence of the air/ground interface 
can be significantly decreased via a data pre-processing procedure. 

Let the ground be 



{(xi,X2) : X2 < a = const.} C K^. 
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Consider a polarized electric field which is generated by a plane wave, initialized at the line 
{x2 — > a, xi £ K} at the moment of time t ^ 0. The following hyperbolic equation can be 
derived from the Maxwell equations in the 2d case 

c{x)utt = Au, {x, t) X (0, oo) . (6.4) 

where the function u{x,t) is a component of the electric field and c{x) := (x) is the spatially 
distributed dielectric constant. We assume that the function c{x) satisfies conditions p.ip in 2d. 
We model imaging of dielectric constants in plastic land mines. In doing so, we do not assume a 
knowledge of the background medium. So, images of land mines are constructed only on the basis 
of values of the dielectric constant c{x) inside of them. 

Let Q be the domain of interest in the ground, where we search for land mines. We set 

n = {{x, y) e (-0.35, 0.35) m x (-0.05, 0.35) m} , 

where "m" stands for meter. Introducing dimensionless spatial variables {x',y') — {x,y) /(0.1m) 
without changing notations for brevity, we obtain the dimensionless domain 

17 = (-3.5,3.5) X (-0.5,3.5). (6.5) 

Hence, the ground is at {x2 = 3.5} and the depth of the domain of interest is 4, which means 40 
cm in real dimensions. Our backreflected signal is measured at the backscattering side, 

the backscattering side is F = {(a;i, X2) : xi e (—0.35, 0.35) , X2 — 3.5} . (6.6) 

It is well known that the maximal depth of an antipersonnel land mine does not exceed about 10 
centimeters. Hence, we model these mines as two small rectangles with the 0.1 m and 0.2 m length 
of sides, and 0.1 m width of sides, respectively. Centers of those rectangles are located at X2 — 2.5 
which is of 10 depth cm in variables with dimensions, see Figure IHTTl 

Tables of dielectric constants [38] show that in the dry sand the dielectric constant Er = 5 and 
Er = 22 in the trinitrotoluene (TNT). Hence, the mine/background contrast is 22/5 « 4. Thus, we 
consider new parameters £j,,t' without changing notations, = er/5,t' = tjy/h. Hence, we obtain 
the following relative values of the dielectric constant in our tests 

c (x) = £r(dry sand) = 1, c(x)= er(TNT) 4. (6.7) 

6.2. Our mathematical model of imaging of explosives hidden in belts worn by 
humans: 3d study. In all places below where the 3d case is discussed, we use the same notation 
for the vector x = (a;, y, z) and for its first coordinate. This does not lead to an ambiguity. In the 
3d case we model the body of a human as a rectangular prism of 2 meters tall, 0.6 meters wide and 
0.16 meters "deep". The vertical coordinate is y and z is responsible for the depth. Hence, in this 
case computational domain G is 

G = {(x, y,z):x£ (-0.5, 0.5) , y G (-1.08, 1.08) , z £ (-0.32, 0.32)} . (6.8) 

We model the belt with explosives as the rectangular prism, which is a subdomain of the first one. 
Sizes of that "belt" are 0.3 meters in the vertical direction, 0.52 meters in horizontal direction and 
0.08 meters of "depth". Hence, dividing by 1 meter, we obtain that these two prisms are respectively 
dimensionless domains and flbeit, 

17 = Gfem = {{x, y,z):x£ (-0.3, 0.3) , y £ (-1, 1) , z e (-0.08, 0.08)} , 17 c G (6.9) 
^beit = {{x,y,z) : X £ (-0. 26, 0.26), 2/ e (-0.15, 0.15) , z e (-0.04,0.04)} C 17. (6.10) 
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On Figure [6?2l the domain G is the largest prism, il is the smaller prism and Albeit is the smallest 
prism. Our backscattering signal is measured at the front side T of the prism 17. The incident plane 
wave propagates along the positive direction of the z— axis. Therefore, the front side of the prism 
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is the backscattering side. We define different boundaries of G and Q. as 



Left side of 51 is = 


{x = -0.3} nH, 


(6.11) 


Right side of £7 is 


{a: = 0.3}nn, 


(6.12) 


Back side of £7 is F;, = 


{z = o.08}nn, 


(6.13) 


Front (backscattering) side of il is F = 


{z = -0.08} nn, 


(6.14) 


Top side of is Ff = 


{y = i}nn, 


(6.15) 


Bottom side of is Thot — 


{y = -i}nn, 


(6.16) 


Front side of G is i9iG = 


{z = -0.32} nG, 


(6.17) 


Back side of G is = 


{z = 0.32}nG, 


(6.18) 


dzG = 


dG\(diG\Jd2G). 


(6.19) 



Therefore, we actually assume here that we measure the backreflected signal at the distance of 
4 cm off the belt. Although this is unrealistic, we can justify this as follows. Suppose that we 
actually measure the backscattering signal on a plane P = {z = — Z, Z > 0.08} . We can approx- 
imately assume that w (cc, y, — 0.08, s) = wq (a^, y, — 0.08, s) for {x^y) G K^\F,s g [s, s]; see sub- 
section 6.5. Recall that the function (f> (x, y, z, s) is the Laplace transform of the data g (x, y, z, t) 
in (|3.4p (subsection 3.2). Using the Green's function for the equation Aw — s^w = in the 
half space {z < —0.08} , we can obtain an integral equation of the first kind with respect to the 
function p (x, y, s) := w (x, y, —0.08, s) . The right hand side of this equation will be the function 
if {x,y, —Z, s) , s will be a parameter and integration will be carried out over the rectangle F. 
This is a convolution equation, which represents a linear ill-posed problem. Algorithms of solving 
convolution equations using the Tikhonov regularization are described in the book [39] • Thus, so- 
lution of this equation would provide us with an approximation of the function w {x, y, —0.08, s) for 
(x, y) G F, s S [s, s] . On the other hand, the latter is the function which we consider as given data 
in our numerical experiments of subsection 6.9. Thus, assuming below that we have the data at F, 
we avoid the intermediate step of solving that integral equation. 

6.3. Data simulation in 2d. To simulate the data for our CIP in 2d, we solve the forward 
problem for equation (|6.4p for the case of the incident plane wave propagating along the negative 
direction of the X2 axis. This plane wave is initialized on the top boundary of the rectangle G of 
Figure 16.11 We simulate the data for the inverse problem using the software package WavES [40] . 
To do that we solve the forward problem via the hybrid FEM/FDM method described in In 
this method the computational domain G is split in two subdomains, G — Gfdm U Gfem, where 

G = [-4,4] X [-1,4],Gfem ■= ^ = (-3.5,3.5) x (-0.5, 3.5) , Gfda/ = G\Gfem, 

see Figure 16.11 Thus the subdomain Gfem := is the same as in (|6.5p . We use structured mesh 
and FDM in Gfdm and non-structured mesh and FEM in Gfem — ^- The space mesh in consists 
of triangles and it consists of squares in Gfdm, with the mesh size h — 0.125 in the overlapping 
regions. At the top and bottom boundaries of G we use first-order absorbing boundary conditions. 
These conditions are exact in our case since we initialize a plane wave in a normal direction to the 
top boundary of G. At the lateral boundaries, the zero Neumann boundary condition is used. Since 
the incident plane wave propagates downwards, then the zero Neumann boundary condition allows 
us to model an infinite space domain in the lateral direction. 
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m>3 


a) t= 3.0 












^1 




b) t= 4.0 




c) t= 5.0 



d) t= 6.0 



Fig. 6.3. Isosurfaces of the computed solution u {x, t) of the forward problem 16.211 . 16.221 for the case of the 
mine-like targets of Figure ISHI for different times. One can observe that values of u {x, t) at the backscattering (top) 
side of the boundary are affected quite significantly by the presence of these targets. Also, values at the bottom 
side are significantly affected. However, since this side is located far away from the backscattering side, then the 
secondary refiected wave does not provide a significant impact on the backscattering side, see d). Values at lateral 
sides are almost the same as ones for the uniform background with c{x) = 1. These observations are the basis for 
our decision about the change of boundary conditions, see 16.25)1 . 



Small square and small rectangle of Figure 16.11 are mine-like targets with c{x) 
them, see ((ST)) . Thus, 



4 inside of 



4 in mine-like targets of Figure 16.1 
1 otherwise. 



(6.20) 



When solving the inverse problem, we assume that the coefficient c{x) is unknown in the rectangle 
C G and has a known constant value c{x) = 1 in G\ft, see Figure [01 The boundary of the 
rectangle G is dG = dGi U 96*2 U dG^. Here, dGi and dG2 are respectively top and bottom sides 
of the largest rectangle of Figure 16.11 and dG^ is the union of left and right sides of this rectangle. 
Let T be the final time for data generation, see the paragraph after (|3.4p in subsection 3.1. We 
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generate the data via solution of the fohowing forward problem 



c{x) utt - Au 
u{x, 0) 



0, ut{x,0) = 0, in G, 
fit), ondGi X (0,ti] 
dtu, on dGi x (ti, T), 
dtu, on dG2 x (0,r), 
0, on dG3 X (0,T), 



0, inGx(0,r), 




(6.21) 



The plane wave with the wave form / (t) is initialized at the top boundary dGi of the compu- 
tational domain G during the time period t e (0, 2t:/lu] = (0, ti], propagates downwards into G and 
is absorbed at the bottom boundary dG2 for all times t e (0,T). In addition, it is also absorbed at 
the top boundary dGi for times t e (27r/a;,T). Here 



We took uj = 7 and T = 6 in (|6.22p for 2d tests. To update tails, we have solved on each iterative 
step the forward problem (|6.2ip . Next, we have calculated the Laplace transform (13.51) to obtain 
the function Wn^i (a;, s) , see p.l3p and (|3.29p . 

The trace g{x,t) of the solution u{x,t) of the forward problem (|6.2ip . (|6.22p is recorded at 
the top boundary T of the domain J7 where we solve the inverse problem, see (|6.6p . This trace 
generates the Dirichlet boundary data ip (x, s) , x ^ T m (I3.16|) (after the Laplace transform). Next, 
the coefficient c{x) is "forgotten", and our goal is to reconstruct this coefficient ior x G ft from the 
data tjj {x, s) . 

6.4. Data simulation in 3d. In this case domains G and f2 are those of (|6.8p and (|6.9p 

respectively. Since the human body consists mostly of water, and the dielectric constant of water 
is about 80 [38j, we set in Test 5 below 



Hence, (|6.23p is a quite heterogeneous and, therefore, a very complicated case. Because of this, we 
start from a simpler problem in our Tests 3,4 via choosing 



As to the subdomain flbeit C fl, we assume that it is filled with an improvised explosive device 
(lED). Analyzing dielectric constants of some materials which might form lEDs |16) . we came to 
the conclusion that we can take c (x) = 3.2 to model an lED. This value of the dielectric constant 
is close to RDX Hexahydro-l,3,5-trinitro-l,3,5-triazine [TB]. Given notations (|6.9p - (|6.19p . (|6.23p . 
(|6.24p . we have simulated the data via solving the forward problem (|6.2ip . (|6.22p . We have used 
the mesh step size h — 0.04 in G. 




(6.22) 




(6.23) 




(6.24) 
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In our 3d tests we took cj = 21 and T = 1 in (|6.22p . To generate backscattered data we solve 
the forward problem (|6.21[) . (|6.22[) in time t = [0, 1] with the time step r = 0.001 using the software 
package WavES [40J. Figure [631 shows isosurfaces of the computed solution u{x,t) of the forward 
problem (|6.2ip . (|6.22p for different times t S (0, 1) for the case when the belt with explosives was as 
the one on Figure [6T2] The trace g {x, t) of the solution u {x, t) of the forward problem ()6.2ip . ()6.22p 
is recorded at the front boundary F of the domain f2, which is the backscattering side of 51, see 
(|6.14p for F. Again, this trace generates the Dirichlet boundary data ijj {x, s) , x G T in p.l6p (after 
the Faplace transform). Next, the coefficient c{x) is "forgotten", and our goal is to reconstruct this 
coefficient for x G f2 from the data ip {x, s) . 

6.5. Boundary conditions on dft\T and the choice of the s interval in 2d. Although 
the above theory requires the knowledge of the function u{x,t) :— g{x,t) at the entire boundary 
(917, the backscattering data are given only on the top part F of the rectangle 51. To see how we 
can complement these data, we analyze the time dependent behavior of the function u (x, t) , which 
is calculated as the solution of the problem ()6.2ip , (|6.22p . Figure 16.31 displays this function for 
different times t € (0, 6) for all a; e fl for the case of two mine-like targets of Figure 16.11 with 
c{x) =4 inside of them and c (a;) = 1 everywhere else, see (|6.7p . We see that values of u (x, t) for 
X £ G are substantially affected by the presence of these inclusions. On the other hand, values at 
lateral sides of the rectangle 51 are affected insignificantly. Values at the lower part of the boundary 
951 are also significantly affected by the presence of those targets. On the other hand, that lower 
part of 951 is located rather far away from the top part of 951. This means that waves reflected 
from the lower part reach the to part of 951 at larger times i > 8. On the other hand, the Laplace 
transform p.Sp actually discounts values of the function u {x, t) for large t: because of the rapid 
decay of the kernel e"'**. 

These observations provide a numerical justification for assigning the following boundary con- 
dition at 951 



Here Wcaic {x, s) is the function w (x, s) which is calculated as the Laplace transform p.Sp of the 
solution of the forward problem (|6.2ip . (|6.22p . On the other hand Wunif {x,s) is the the Laplace 
transform of the solution of this problem for the uniform medium with c (x) = 1. 

Consider now the function , a; G 51. Figure [^13] displays graphs of the function q (x, s) along the 
top boundary of the rectangle 51 for different values of the pseudo frequency s. One can observe 
that for s = 2,3,4,5,7 each graph has dents. The locations of these dents exactly correspond to 
projections of two targets of Figure [Ql on the top boundary of 51. Therefore, values of the function 
q{x, s) ,x € F carry an information about the horizontal coordinate of this inclusion. However, 
figuring out vertical coordinates of targets is a more difficult task. To do this, one needs to apply 
the above algorithm. We also observe that values of \q {x, s)\ for s > 7 are much lower than those 
for s < 5. Therefore, to solve the inverse problem in 2d, we have chosen the s— interval as 



6.6. 3d case: boundary conditions on 951\F and the choice of the s— intervaL Just 
as in the 2d case, we have chosen boundary conditions as in (j6.25p . To justify this, we present 
Figures [6.6116.91 One can see from these figures that values of the function w{x,s) at all parts of 
951\F, except of the back side F;, of the prism 51, are about the same as ones for the case c (x) = 1, 




(6.25) 



s e [2,3] ,/i = 0.05. 



(6.26) 
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which is the value of this coefficient outside of our domain of interest ft. As to the surface Fb, it 
corresponds to the transmitted signal and values of w (x, s) here are far from those of the uniform 
background outside of 17. However, just as in the 2d case, the transmitted side Ti, is located far 
from the backscattering side F. Therefore, the Laplace transform p.5|) diminishes the influence of 
waves reflected from Ff,, at least for large values of the parameter s, see Figures 16^16.71 Hence, the 
amplitude of reflected waves from this side is small when they reach F, compared with reflections 
at F from the target ilbeit- This provides a numerical justiflcation of (|6.25p in the 3d case. We call 
the resulting boundary function -0 (x, s) "immersed boundary data", see Figure [6. 11 1 

6.7. Numerical tests for the 2d case. Let Ucaic {x,t) , x e F be the calculated solution 
of the forward problem (|6.2ip . (|6.22l) at the backscattering side F of the boundary dft. We have 
introduced a random noise in the function Ucaic {x,t) ,x Cz T as 

= u,aic [1 + a^a (u^,^ - , (6.27) 

where (x*-*-* , t*-^-*) G F x (0,T) are mesh points, itmax and are maximal and minimal values of 
Ucaic {x, t) for cc G F, numbers aj G (—1, 1) are randomly distributed and a = 0.05. Thus, the noise 
level was 5%. In both Test 1 and Test 2 the correct coefficient c{x) is the same as in (|6.20p . 

Test 1. In this test the initial guess for the tail function Vi.i the function was computed via 
p.l3p for the case of the homogeneous domain G with c{x) = 1. Next, the algorithm of subsection 
3.4 was applied to reconstruct the true function c{x) in (j6.20p . The computed image is presented 
on Figure l6.12l -a). We observe that both the location and the contrast of both mine-like targets 
are reconstructed accurately. The number of inner iterations with respect to tails was m — 5. The 
stopping criterion ()6.2p - (|6.3p was achieved at 03,5 {x) , i.e. we have stopped at n = 3. By (|6.26p . 
(I3.18P and (|3.19p this corresponds to s G [2.85, 2.90] . The reconstructed dielectric constant in this 
test is C3^5 (x) = 4.07 inside of both imaged mine-like targets and 03^5 (x) = 1 at all other points 
of J7. To see what happens in an ideal case when the exact tail function V* {x) is known, we 
refer to Figure [6.121 b). which corresponds to the function ci 1 (x) . Figure [6.121 b) confirms that the 
reconstruction is perfect in this case. 

Test 2. In this test we choose the initial guess for the tail function Vi^ias an initial guess for 
the tail function Vi^i we take the function computed via p.4ip and use the algorithm of subsection 
3.4 to reconstruct the dielectric constant of Figure [6.131 a). In this test the reconstructed dielectric 
constant is cg^e (x) = 4.27 inside mine-like targets and cg^e (x) = 1 at all other points of fl. This 
reconstruction was obtained on the pseudo-frequency interval s G [2.6, 2.65] and after 6 iterations 
with respect to the tail function. In other words, we took the number of inner iterations with 
respect to tails m = 6 and the stopping criterion (|6.2p - (|6.3p was achieved at n — 8. 

Thus, in both Tests 1,2 reconstructions were accurate ones. 

6.8. Numerical test for the 3d case. In Test 3 and Test 4 we present results for the case 
(|6.24p . and in Test 5 - for the case ()6.23p . The same random noise of 5% was introduced as the one 
in (jOTi) . 

Test 3. In this test we took the first guess for the tail function Vi.i (x) the same as the one 
for the uniform background when c (x) = 1 for x £ G. Using Figures 16.61 16.71 and analyzing the 
backscattered data 'ip{x^ s) for x G F we have decided to choose the interval of pseudo frequencies 
as 



sG [4, ll],/i = l. 
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Then we have used the algorithm of subsection 3.4 to reconstruct the dielectric constant in the 
belt of Figure 16.21 -c) . Figure 16.141 -al presents reconstruction of the dielectric constant c{x) for the 
unrealistic case when we know the exact tail function. In this case we observe that the reconstruction 
is perfect. 

Figure 16.141 -b) shows the maximal values of the reconstructed function c{x) when the initial 
tail Vi^i (x) was computed from the homogeneous domain with c{x) = 1 for x ^ G. We observe 
that the location and the contrast of the explosive-like target are reconstructed accurately. The 
reconstructed dielectric constant in this test is 03^2 (x) = 3.54 inside the "belt", and c{x) ~ 1 at 
all other points of il. We took the number of inner iterations with respect to tails m — 2, and the 
stopping criterion (|6.2p - (|6.3p was achieved at n = 3, which corresponds to s S [8,9] — [54,53] in 
p.lSp . (I3.19p . We conclude that this reconstruction is accurate. 

Test 4. In this test we took the tail Vi^i (x) the same as in our above theory, see (|3.39p . p.40p . 
p.4ip . Analyzing results of Test 3 we have also decided to refine the pseudo- frequency interval in 
this test. Indeed, we got our final image of Test 3 for s G [8, 9] . Hence, we decided to take the 
interval of pseudo frequencies 

s € [8,8.85] ,/i = 0.05. (6.28) 

Next, we have used the algorithm of subsection 3.4. We took 171 — 2. The stopping criterion 
()6.2p - (|6.3p was achieved at n = 3, which corresponds to s G [8.7,8.75] — [34,33] in p.l8p . (|3.19p . 
The reconstructed function 03^2(2;) is depicted on Figure [6.151 -a). In this test the reconstructed 
dielectric constant is 03^2 (x) = 3.54 inside the belt and C3.2 (x) = 1 at all other points of il. Thus, 
the reconstruction was again a quite accurate one. 

Test 5. This is the most challenging test, because the medium is a quite heterogeneous one: 
there are substantial contrasts between three values of the target function c{x). Indeed, we have used 
the case (this test we took the first guess for the tail Vi^i (x) as in our above theory see (|3.39p . 
f (|3.4ip . i.e. the same as the one in Test 4. Used results of Test 4 we took now the interval of 
pseudo- frequencies 

3 £ [8.0,8.8] ,/i = 0.05. 

Next, we have used the algorithm of subsection 3.4. However, since we know that the dielectric 
constant of the human body is large, c — 80, then we have truncated to 1 those values of computed 
functions c„^i (x) , which exceeded 10. In other words ()6.ip . was replaced with 

_ f (^n,i (x) if Cn,i (x) e [1, 10] and x e Q, _ 

1^ 1 if either Cn,i (x) < 1, or Cn,i (x) > 10, or x ^ fl. 

We took m = 3 and the stopping criterion (|6.2p - (l6.3p was achieved at n = 2. The latter corresponds 
to s G [8.70,8.75] = [53,52] in (|3.18p . (|3.19p . The reconstructed function 02,3(2:) is depicted on 
Figure [B.15l -b]l. The reconstructed dielectric constant is c = 3.09 inside the belt, and c = 1 at all 
other points of Q. Therefore, the reconstruction is again a quite accurate one even in this most 
difficult case. 

7. Summary. We have presented a new approximate mathematical model. This model 
amounts to the truncation of the asymptotic series with respect to 1/s, where s >> 1 is the 
upper limit of the positive parameter of the Laplace transform of the solution of the Cauchy prob- 
lem p.2p . (|3.3p . However, this truncation is done only on the first iteration of our method to ensure 
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estimate (|5.36p for the accuracy of the first tail function VVi i . No other "special" simplifying 

l+a 

assumptions are made. On the basis of this new model, we have developed a new convergence 
analysis, which is more realistic than the one of our first publications [21 [8] about this topic. This 
time we estimate tail functions. Tails were not estimated in our previous publications, and this is 
a significantly new element of the convergence analysis here. 

We have modified our approximately globally convergent algorithm for the case of backscatter- 
ing data. To do so, we have used a computational observation that one can replace the unknown 
Dirichlet boundary condition on the non-backscattering part of the boundary with the data ob- 
tained for the case of the uniform background, which is assumed to be known outside of the domain 
of interest (but not inside of it), see (|6.25|) . Therefore, our previously developed technique for the 
case when the Dirichlet data are given at the entire boundary, works. Our numerical tests confirm 
this. 

Our numerical tests 1-4 demonstrate that the case when the first tail is taken the same as the 
one for the uniform medium with c (x) = 1 provides almost the same results as ones for the new 
tail function. Numerical studies demonstrate the accuracy of our technique. It is worthy to note 
that we have obtained an accurate image even in the most difficult case of Test 5 when the medium 
was quite a heterogeneous one, see ()6.23|) . 

We believe that results of Tests 1-5 combined with results for blind experimental data of |29[ [5U] 
and section 6.9 of [6] confirm the validity of our approximate mathematical model, as indicated in 
Steps 4-6 of section 2. 
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a) s = 7 b) s = 19 

Fig. 6.4. Backscattered data at the top boundary T of the function q(x, s) at different values pseudo-frequencies 
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a) h) t = 0.85 c) t = 0.95 



Fig. 6.5. Isosurfaces of the computed solution u{x,t) of the forward problem (|6.2ip . (|6.22p at different 
times t with the plane wave initialized at the front boundary of G on the mesh with the mesh size h — 0.04. 
Test was performed in time t — [0, 1] with time step r = 0.001. 





a) s = 8 b) s = 7 

Fig. 6.6. Backscattered data ^p{x,s),x G F at different pseudo-frequencies s £ [7; 20]. 
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Fig. 6.7. Backscattered data ^p{x,s),x £ F at diflFerent pseudo-frequencies s £ [1,6]. One can see from 
Figures [6.61 and 1 6.71 that one should take s > 4 in the reconstruction algorithm in 3d tests. 
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c) 

Fig. 6.8. Analysis of the scattered data at the left F; and right Tr boundaries of the Gfeai domain: 
a) scattered data %l){x,s),x £ Fj superimposed with the scattered data tp{x,s),x £ Vr- one can see that 
the dent is very small. The data are scattered from the belt modeling and explosive, see Figure [62] Here, 
c = 3.2 inside the belt and c = 1 at all other points of Gfem', b) homogeneous data s) with c = 1 in 
Gfeai; c) Superimposed data of a) and b). 
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1 



b) 




c) 

Fig. 6.9. Analysis of the scattered data at the top Ft and bottom Tbot boundaries of the Gfem domain: 
a) scattered data ip(x, s),x £Vt superimposed with scattered data ip{x, s),x £ Ftot- The data are scattered 
from the belt with explosive of Figure 16.21 with c = 3.2 inside the belt and c = 1 at all other points of 
Gfem; b) homogeneous data xp{x,s) with c = 1 in Gfem', c) Superimposed data of a) and b). 
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Fig. 6.10. Test 5. Analysis of scattered data at the top Ft and bottom Ti,ot boundaries of the Gfem' 
a) scattered data tp{x, s),x £Tt superimposed with scattered data ^p{x, s), x £ Ftot- b) Superimposed data 
of a) (low figure) and homogeneous data tp{x,s) (top figure) computed with c = 1 in Gfem- 




Fig. 6.11. Test 5. Backscattered data ^l){x, s),x £T immersed into data ^p{x, s), x € dQ,\T computed 
with c = 1 in f2 = Gfem- Data are presented at different pseudo-frequencies s = 6, 7, 8, 9. 
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Fig. 6.12. Test 1: a) The computed function ca.s (x) . Maximal values of this function are 4.07 in both imaged 
mine-like targets and 03,5 {x) = 1 outside of imaged targets. The image is accurate: compare with Figure [6TT}c) and 
writh 1 16.201 . b) The reconstructed function ci 1 (x) for the the case when the exact tail function V {x) is known. 
However, this is an unrealistic case, which is presented here only to demonstrate the accuracy of our method in the 
ideal case. 




a) 



b) 



Fig. 6.13. Test 2: a) The computed function cg g {^) ■ Maximal values of this function are 4.27 in both imaged 
mine-like targets and cg^g (x) = 1 outside of imaged targets. The image is accurate: compare with Figure IHT^ c') and 
with II6.2OI I. b) The reconstructed function ci.i {x) for the the case when the exact tail function V (x) is known. 



44 



L.BEILINA AND M.V.KLIBANOV 




a) maxci.i « 3.2 




b) niaxc3,2 ~ 3.54 

Fig. 6.14. Test 3. Reconstruction of a belt with explosives for the case (|6.24ll . a) The computed 
function ci,i (a;) for the case when the exact tail function V* (x) is known. The reconstruction is perfect. 
However, this is an unrealistic scenario. We display it here only to show that our method is accurate in 
an ideal case, b) The reconstruction for the case when the initial tail function Vi,i (x) is taken the same 
as the one for the uniform medium with c{x) = 1. The computed function £3,2 (x) is shown. Observe that 
maxc3,2 (x) = 3.54 inside of the imaged inclusion (belt with explosives) and 03,2 (a:) = 1 everywhere else. 
The image is accurate. 
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b) maxc2,3 ~ 3.09 

Fig. 6.15. Reconstruction with back-scattered data, a) Test 4. The computed image of the function 
C3,2. Observe that max £3,2 (x) = 3.54 inside of the imaged inclusion (belt with explosives) and 03,2 (x) = 1 
everywhere else, b) Test 5. This is the most challenging test since the medium is a quite heterogeneous 
one, see (|6.23p . The computed image of the function £2,3 (x) . One can see that max £2,3 (x) — 3.09 inside 
of the imaged inclusion ("belt" with explosives) and C2,3{x) — 1 everywhere else. Both images are quite 
accurate ones. 



